# -------------------------------------------------------------------------
#
# Script to reproduce results in 
#
  # "The role of grandiose and vulnerable narcissism on mask 
  #   wearing and vaccination during the COVID‑19 pandemic" 
#
# Journal: Current Psychology
#
# Authors: Z Fazekas & P Hatemi
#
# Content: main paper (with some additional material reported in SI)
#
# Please find sessionInfo() at the end of this script for last run instance
#
# -------------------------------------------------------------------------

# packages
# note: install before, if needed
library("tidyverse")
library("psych")
library("lavaan")
library("Hmisc")
library("texreg")
library("broom")
library("survey")

library("conflicted")
conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")
conflict_prefer("summarize", "dplyr")

two_sd <- function(x){
  (x - mean(x, na.rm = TRUE))/(2*sd(x, na.rm = TRUE))
}

us <- read_csv("cp-yougov_narc.csv")

# NPI (25 items) ----------------------------------------------------------
npi_vars   <- paste0("PI", 1:25)
# narcissistic answers as top
to_reverse <- c(1, 6, 7, 8, 9, 13, 14, 
                16, 18, 19, 21, 22, 23, 24)  # from NPI
us[, npi_vars] <- us[, npi_vars] - 1
us[, npi_vars[to_reverse]] <- 1 - us[, npi_vars[to_reverse]] 

# relabel to original npi item number
fullnpi_vars <- paste0("npi_",
                       c(1, 4, 5, 7, 10, 11, 12, 13, 14, 15, 19, 20, 
                         24, 25, 26, 27, 28, 29, 30, 32, 33, 34, 36, 38, 40))

# use different rename
names(us)[which(names(us) %in% npi_vars)] <- fullnpi_vars

## 3-factor variables
lead_auth   <- c("npi_1", "npi_5", "npi_10", "npi_11", 
                 "npi_12", "npi_27", "npi_32", "npi_33", "npi_34", "npi_36", "npi_40")
grand_exhib <- c("npi_4", "npi_7", "npi_15", "npi_19", 
                 "npi_20", "npi_26", "npi_28", "npi_29", "npi_30", "npi_38")
ent_exp     <- c("npi_13", "npi_14", "npi_24", "npi_25")
npi25_vars <- c(lead_auth, grand_exhib, ent_exp)

psych::alpha(us[, lead_auth])   ## 0.74 reliability
psych::alpha(us[, grand_exhib]) ## 0.75 reliability
psych::alpha(us[, ent_exp])     ## 0.48 reliability
psych::alpha(us[, npi25_vars])  ## 0.82 reliability

## 3-factor CFA check
## reported in main paper
## ((DWLS estimator, robust: 0.899 (CFI), 0.889 (TFI), 
## RMSEA of 0.049 (90% CI, 0.046, 0.052), and SRMR 0.09).)
npi.3factor <-  'lead_auth =~ npi_1 + npi_5 + npi_10 + npi_11 + npi_12 + 
                              npi_27 + npi_32 + npi_33 + npi_34 + 
                              npi_36 + npi_40
                 grand_exh =~ npi_4 + npi_7 + npi_15 + npi_19 + npi_20 + 
                              npi_26 + npi_28 + npi_29 + npi_30 + npi_38
                 ent_exp   =~ npi_13 + npi_14 + npi_24 + npi_25'
npi.3factor.fit <- cfa(npi.3factor, data = us, ordered = npi25_vars)
summary(npi.3factor.fit, fit.measures = TRUE, standardized = TRUE)

# create variables
us$npi25      <- rowSums(us[, npi25_vars])/length(npi25_vars)
us$leadauth   <- rowSums(us[, lead_auth])/length(lead_auth)
us$entexp     <- rowSums(us[, ent_exp])/length(ent_exp)
us$grandexhib <- rowSums(us[, grand_exhib])/length(grand_exhib)

# store out narcissism variable names
narc25_vars <- c("npi25", "leadauth", "grandexhib", "entexp")

# HSNS --------------------------------------------------------------------
hsns_orig <- paste0("PI2_q", 26:35)
hsns_vars <- paste0("hsns_", 1:10)

names(us)[which(names(us) %in% hsns_orig)] <-  hsns_vars
us[, hsns_vars] <- apply(us[, hsns_vars], 2, function(x) {(x - 1)/4})

# using Fossati et al model
## 2-factor CFA
## 0.901 (CFI), 0.868 (TFI), 
## RMSEA of 0.066 (90% CI, 0.058, 0.075), and SMRS 0.050.
hsns.2factor <-   'ego_cent =~ hsns_4 + hsns_5 + hsns_8 + hsns_10
                 judge_over =~ hsns_1 + hsns_2 + hsns_3 + hsns_6 + 
                               hsns_7 + hsns_9'
hsns.2factor.fit <- cfa(hsns.2factor, data = us)
summary(hsns.2factor.fit, fit.measures = TRUE, standardized = TRUE)

# create variables
ego_cent     <- c("hsns_4", "hsns_5", "hsns_8", "hsns_10")
judge_over   <- c("hsns_1", "hsns_2", "hsns_3", "hsns_6", "hsns_7", "hsns_9")
us$egocent   <- rowSums(us[, ego_cent])/length(ego_cent)
us$judgeover <- rowSums(us[, judge_over])/length(judge_over)
us$hsns_all <- rowSums(us[, hsns_vars])/length(hsns_vars)

psych::alpha(us[, hsns_vars]) # 0.73
psych::alpha(us[, ego_cent])   # 0.63 
psych::alpha(us[, judge_over]) # 0.69
hsns_varnames <- c("hsns_all", "egocent", "judgeover") 

all_narc <- c(narc25_vars, hsns_varnames)

apply(us[, all_narc], 2, function (x) {mean(x, na.rm = TRUE)})
apply(us[, all_narc], 2, function (x) {sd(x, na.rm = TRUE)})

# Table 3 (unweighted)
round(cor(us[, all_narc], use = "complete"), 3)[1:4, 5:7]
round(rcorr(as.matrix(us[, all_narc]))$P[1:4, 5:7], 3)

# Recode covariates -------------------------------------------------------
# states with mask mandate added
us <- us |> 
  mutate(dem_rep = pid7 - 1,
         dem_rep = ifelse(pid7 == 8, 3, dem_rep),
         female  = gender - 1,
         age     = 2021 - birthyr,
         age     = age - min(age),
         not_caucasian = ifelse(race == 1, 0, 1),
         edu     = educ - 1,
         state_mask = ifelse(inputstate == 2  |
                               inputstate == 4  |
                               inputstate == 12 |
                               inputstate == 13 |
                               inputstate == 16 |
                               inputstate == 19 |
                               inputstate == 28 |
                               inputstate == 29 |
                               inputstate == 30 |
                               inputstate == 31 |
                               inputstate == 38 |
                               inputstate == 40 |
                               inputstate == 45 |
                               inputstate == 46 |
                               inputstate == 37 |
                               inputstate == 48 |
                               inputstate == 56, 0, 1))

covs <- c("dem_rep", "edu", "female", "age", "not_caucasian", "state_mask")
round(rcorr(as.matrix(us[, c(covs, all_narc)]))$r[1:6, 7:13], 3)

# COVID worry (predictor) -------------------------------------------------
us$covid_worry       <- 5 - us$Q38
round(rcorr(as.matrix(us[, c("covid_worry", all_narc)]))$r[1, 2:8], 3)
round(rcorr(as.matrix(us[, c("covid_worry", all_narc)]))$P[1, 2:8], 3)

all_pred <- c(covs, "covid_worry")
# Outcomes of interest ----------------------------------------------------
us <- us |> 
  mutate(covid_mask_own    = 2 - Q39,
         covid_mask_others = 2 - Q40,
         covid_mask_told   = 2 - Q41, # NA is for those who do not wear as well
         covid_vaccine_2 = case_when(Q42 == 2 ~ 0, # no
                                   Q42 == 4 ~ NA_real_, # dk
                                   Q42 == 1 | Q42 == 3 ~ 1 # yes
         ))
covid_vars <- names(us)[str_detect(names(us), "covid_")]
covid_vars <- covid_vars[str_detect(covid_vars, "worry", 
                                    negate = TRUE)]

# final selection of data (for regressions)
us <- select(us, c("caseid", "weight", 
                   all_of(covid_vars), 
                   all_of(all_pred),
                   all_of(all_narc)))

# cor.test(us$covid_worry, us$npi25, use = "complete")
# cor.test(us$covid_worry, us$leadauth, use = "complete")
# cor.test(us$covid_worry, us$entexp, use = "complete")
# cor.test(us$covid_worry, us$grandexhib, use = "complete")
# cor.test(us$covid_worry, us$hsns_all, use = "complete")
# cor.test(us$covid_worry, us$egocent, use = "complete")
# cor.test(us$covid_worry, us$judgeover, use = "complete")

# distribution plots
grand_narc <- all_narc[1:4]
vuln_narc  <- all_narc[5:7]

us |> 
  select(all_of(grand_narc)) |> 
  pivot_longer(cols = all_of(grand_narc)) |> 
  mutate(name_s = case_when(name == "npi25" ~ "1. Full narcissism (NPI)", 
                            name == "leadauth" ~ "2. Leadership/Authority",
                            name == "grandexhib" ~ "3. Grandiose exhibitionism",
                            name == "entexp" ~ "4. Entitlement/Exploitativeness")) |> 
  ggplot(aes(x = value)) +
  geom_histogram() +
  facet_wrap(~name_s, scale = "free_y") +
  theme_minimal() +
  labs(y = "Count", x = "") +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom")

# uncomment to save
# ggsave(file = "distr-grand.pdf",
#        width = 12, height = 6,
#        device = cairo_pdf)


us |> 
  select(all_of(vuln_narc)) |> 
  pivot_longer(cols = all_of(vuln_narc)) |> 
  mutate(name_s = case_when(name == "hsns_all" ~ "1. HSNS full",
                            name == "egocent" ~ "2. Egocentrism",
                            name == "judgeover" ~ "3. Oversensitivity\nto judgement")) |>
  ggplot(aes(x = value)) +
  geom_histogram() +
  facet_wrap(~name_s, scale = "free_y") +
  theme_minimal() +
  labs(y = "Count", x = "") +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom")

# uncomment to save
# ggsave(file = "distr-vuln.pdf",
#        width = 12, height = 6,
#        device = cairo_pdf)


# keeping only essential
rm(list = setdiff(ls(), 
                  c("us", "covid_vars", "all_pred", "all_narc", "two_sd")))

# recode continuous for easier comparisons
us_reg <- us |> 
  mutate_at(c("dem_rep", "edu", "age", "covid_worry",
              all_narc), two_sd)

# create weighted design datasets for later
us_weighted <- svydesign(ids = ~1, data = us, weights = ~weight)
us_reg_weighted <- svydesign(ids = ~1, data = us_reg, weights = ~weight)

# table 1 (weighted)
apply(us[, covid_vars], 2, 
      function (x){round(svymean(x, design = us_weighted, na.rm = TRUE), 3)})
# table(us[, covid_vars[2]]) # 900 said yes/1

# table 2 (weighted) + unweighted alpha calculated above (variable grouping dropped)
apply(us[, all_narc], 2, 
      function (x){round(svymean(x, design = us_weighted, na.rm = TRUE), 2)})
apply(us[, all_narc], 2, 
      function (x){round(sd(x, na.rm = TRUE), 2)})

# in text predictor weighted averages
# note: age is already rescaled here, for original run prior check
apply(us[, all_pred], 2, 
      function (x){svymean(x, design = us_weighted, na.rm = TRUE)}) 
apply(us[, all_pred], 2, 
      function (x){sd(x, na.rm = TRUE)})
# Note: covid_worry correlations for footnote calculated before

# All unweighted regression and figure -----------------------------------
# main text
us_reg <- us_reg |> 
  select(all_of(covid_vars), all_of(all_pred), all_of(all_narc)) |> 
  pivot_longer(cols = all_of(covid_vars),
               names_to = "outcomes", values_to = "y")

us_reg <- us_reg  |>  
  mutate(outcomes = ifelse(outcomes == "covid_mask_others", "covid_mask_rest", outcomes))

models <- us_reg |>   
  nest_by(outcomes) |>
  mutate(mod_npi    = list(glm(y ~ npi25, data = data,
                               family = binomial(link = "logit"))),
         mod_npi_wc = list(glm(y ~ npi25 + covid_worry + dem_rep + 
                                edu + female + age + not_caucasian + 
                                factor(state_mask), data = data,
                               family = binomial(link = "logit"))),
         mod_npif    = list(glm(y ~ leadauth + grandexhib + entexp, data = data,
                                family = binomial(link = "logit"))),
         mod_npif_wc = list(glm(y ~ leadauth + grandexhib + entexp + 
                                 covid_worry + dem_rep +
                                edu + female + age + not_caucasian + 
                                factor(state_mask), data = data,
                                family = binomial(link = "logit"))),
         mod_hsns    = list(glm(y ~ hsns_all, data = data,
                                family = binomial(link = "logit"))),
         mod_hsns_wc = list(glm(y ~ hsns_all + covid_worry + dem_rep +
                                 edu + female + age + not_caucasian + 
                                 factor(state_mask), data = data,
                                family = binomial(link = "logit"))),
         mod_hsnsf    = list(glm(y ~ egocent + judgeover, data = data,
                                 family = binomial(link = "logit"))),
         mod_hsnsf_wc = list(glm(y ~ egocent + judgeover + covid_worry + dem_rep +
                                 edu + female + age + not_caucasian + 
                                 factor(state_mask), data = data,
                                 family = binomial(link = "logit"))))

reg_plot <- bind_rows(data.frame(models |> summarise(broom::tidy(mod_npi)),
                     model = "npi25"),
          data.frame(models |> summarise(broom::tidy(mod_npi_wc)),
                     model = "npi25_wc"),
          data.frame(models |> summarise(broom::tidy(mod_npif)),
                     model = "npi25f"),
          data.frame(models |> summarise(broom::tidy(mod_npif_wc)),
                     model = "npi25f_wc"),
          data.frame(models |> summarise(broom::tidy(mod_hsns)),
                     model = "hsns25"),
          data.frame(models |> summarise(broom::tidy(mod_hsns_wc)),
                     model = "hsns25_wc"),
          data.frame(models |> summarise(broom::tidy(mod_hsnsf)),
                     model = "hsns25f"),
          data.frame(models |> summarise(broom::tidy(mod_hsnsf_wc)),
                     model = "hsns25f_wc"))

res_unweighted <- reg_plot |> 
  filter(term %in% all_narc) |> 
  select(outcomes, term, estimate, std.error, model, p.value) |> 
  mutate(wc = ifelse(str_detect(model, "wc"), "with controls", "without controls"),
         term_s = case_when(term == "npi25" ~ "NPI (full)",
                            term == "leadauth" ~ "Leadership-\nAuthority",
                            term == "grandexhib" ~ "Grandiose\nExhibitionism",
                            term == "entexp" ~ "Entitlement-\nExploitativeness",
                            term == "hsns_all" ~ "HSNS (full)",
                            term == "egocent" ~ "Egocentrism",
                            term == "judgeover" ~ "Judgement\noversensitivity"),
         term_f = factor(term_s,
                         levels = c("Judgement\noversensitivity",
                         "Egocentrism",
                         "HSNS (full)",
                         "Entitlement-\nExploitativeness",
                         "Grandiose\nExhibitionism",
                         "Leadership-\nAuthority",
                         "NPI (full)")),
         outcomes_s = case_when(outcomes == "covid_mask_own" ~ "(1) Wears mask",
                  outcomes == "covid_mask_rest" ~ "(2) Others should wear mask",
                  outcomes == "covid_mask_told" ~ "(3) Told others to wear mask",
                  outcomes == "covid_vaccine_2" ~ "(4) Get vaccinated"),
         is_sig = ifelse((estimate + 1.96*std.error) * (estimate - 1.96*std.error) > 0 ,
                         1, 0),
         p_stars = case_when(p.value <= 0.05 & p.value > 0.01 ~ "*",
                             p.value <= 0.01 & p.value > 0.001 ~ "**",
                             p.value <= 0.001 ~ "***"),
         p_stars = ifelse(is.na(p_stars), "", p_stars))

ggplot(res_unweighted, aes(x = term_f, y = estimate,
             ymin = estimate - 1.96*std.error,
             ymax = estimate + 1.96*std.error, shape = wc,
             colour = factor(is_sig),
             label = p_stars)) +
  geom_hline(yintercept = 0, colour = "grey80", alpha = 0.5, linetype = 2) +
  scale_colour_manual(values = c("grey80", "black"), guide = "none") +
  geom_linerange(position = position_dodge(width = 0.35)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  facet_wrap(~outcomes_s, ncol = 4) + 
  geom_text(aes(y = -1.25), position = position_dodge(width = 0.35)) +
  coord_flip() +
  labs(x = "", y = "Effect associated with 2SD difference (logit)",
       shape = "Model", caption = "Whiskers are 95% confidence intervals.\nStatistically significant (p < 0.05) coefficients in black.\n*** p < 0.001, ** p < 0.01, * p < 0.05.") +
  theme_minimal() +
  scale_shape_manual(values = c(16, 15)) +
  theme(legend.position = c(0.48, 0.9)) +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom",
        legend.direction = "horizontal")

## uncomment to save
# ggsave(file = "model-coefs.pdf",
#        width = 12, height = 6,
#        device = cairo_pdf)

## uncomment to save
# export the regression tables (html format)
# htmlreg(models$mod_npi,
#           caption = "",
#           custom.model.names = c("Wears mask", "Others wear",
#                                  "Told others", "Get vaccinated"),
#           include.deviance = FALSE,
#           include.bic  = FALSE,
#           custom.coef.names = c("Intercept", "Narcissism: NPI25 (full)"),
#           custom.gof.names = c("AIC", "LogLik", "N"), 
#           file = "npi25-models.html")
# htmlreg(models$mod_hsns,
#           caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#           include.bic  = FALSE,
#           custom.coef.names = c("Intercept", "Narcissism: HSNS (full)"),
#           custom.gof.names = c("AIC", "LogLik", "N"), 
#           file = "hsns-models.html")
# htmlreg(models$mod_npif,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: Leadership-authority",
#                               "Narcissism: Grandiose exhibitionism",
#                               "Narcissism: Entitlement-exploitativeness"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "npi25-sub-models.html")
# htmlreg(models$mod_hsnsf,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: Egocentrism", 
#                               "Narcissism: Oversensitivity to judgement"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "hsns-sub-models.html")
# 
# htmlreg(models$mod_npi_wc,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: NPI25 (full)",
#                               "Worried about COVID-19",
#                               "Party ID (Dem to Rep)",
#                               "Education",
#                               "Sex (Female = 1)",
#                               "Age",
#                               "Not caucasian",
#                               "Mask mandate in state"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "npi25-models-wc.html")
# htmlreg(models$mod_hsns_wc,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: HSNS (full)",
#                               "Worried about COVID-19",
#                               "Party ID (Dem to Rep)",
#                               "Education",
#                               "Sex (Female = 1)",
#                               "Age",
#                               "Not caucasian",
#                               "Mask mandate in state"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "hsns-models-wc.html")
# htmlreg(models$mod_npif_wc,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: Leadership-authority",
#                               "Narcissism: Grandiose exhibitionism",
#                               "Narcissism: Entitlement-exploitativeness",
#                               "Worried about COVID-19",
#                               "Party ID (Dem to Rep)",
#                               "Education",
#                               "Sex (Female = 1)",
#                               "Age",
#                               "Not caucasian",
#                               "Mask mandate in state"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "npi25-sub-models-wc.html")
# htmlreg(models$mod_hsnsf_wc,
#         caption = "",
#         custom.model.names = c("Wears mask", "Others wear",
#                                "Told others", "Get vaccinated"),
#         include.deviance = FALSE,
#         include.bic  = FALSE,
#         custom.coef.names = c("Intercept", "Narcissism: Egocentrism", 
#                               "Narcissism: Oversensitivity to judgement",
#                               "Worried about COVID-19",
#                               "Party ID (Dem to Rep)",
#                               "Education",
#                               "Sex (Female = 1)",
#                               "Age",
#                               "Not caucasian",
#                               "Mask mandate in state"),
#         custom.gof.names = c("AIC", "LogLik", "N"), 
#         file = "hsns-sub-models-wc.html")


# for reusing script
all_pred <- c(all_pred[7], all_pred[1:6])

## Run and export the survey weights model
mods <- data.frame(outcomes = c(rep(covid_vars, 4), rep(covid_vars, 4)),
                   model = c(rep("without controls", 16), rep("with controls", 
                                                              16)), 
                   preds = c(rep("npi", 4), rep("npi_sub", 4), rep("hsns", 4), 
                             rep("hsns_sub", 4),
                             rep("npi", 4), rep("npi_sub", 4), rep("hsns", 4), 
                             rep("hsns_sub", 4)),
                   form  = c(paste0(covid_vars, " ~ ", all_narc[1]),
                             paste0(covid_vars, " ~ ", paste(all_narc[2:4], 
                                                             collapse = " + ")),
                             paste0(covid_vars, " ~ ", all_narc[5]),
                             paste0(covid_vars, " ~ ", paste(all_narc[6:7], 
                                                             collapse = " + ")),
                             paste0(covid_vars, " ~ ", all_narc[1], 
                                    " + ", paste(all_pred, collapse = " + ")),
                             paste0(covid_vars, " ~ ", paste(all_narc[2:4], 
                                                             collapse = " + "), 
                                    " + ", paste(all_pred, collapse = " + ")),
                             paste0(covid_vars, " ~ ", all_narc[5], 
                                    " + ", paste(all_pred, collapse = " + ")),
                             paste0(covid_vars, " ~ ", paste(all_narc[6:7], 
                                                             collapse = " + "), 
                                    " + ", paste(all_pred, collapse = " + "))))

results <- data.frame()
all_mods <- list()
for (i in 1:nrow(mods)){
  m <- svyglm(formula(mods$form[i]), design = us_reg_weighted, 
              family = quasibinomial(link = "logit"))
  all_mods[[i]] <- m
  temp <- tidy(m)
  temp[, names(mods)] <- mods[i, ]
  results <- rbind(results, temp)
}
results <- select(results, -form)

res_weighted <- results |> 
  filter(term %in% all_narc) |> 
  select(outcomes, term, estimate, std.error, model) |> 
  mutate(wc = model,
         term_s = case_when(term == "npi25" ~ "NPI (full)",
                            term == "leadauth" ~ "Leadership-\nAuthority",
                            term == "grandexhib" ~ "Grandiose\nExhibitionism",
                            term == "entexp" ~ "Entitlement-\nExploitativeness",
                            term == "hsns_all" ~ "HSNS (full)",
                            term == "egocent" ~ "Egocentrism",
                            term == "judgeover" ~ "Judgement\noversensitivity"),
         term_f = factor(term_s,
                         levels = c("Judgement\noversensitivity",
                                    "Egocentrism",
                                    "HSNS (full)",
                                    "Entitlement-\nExploitativeness",
                                    "Grandiose\nExhibitionism",
                                    "Leadership-\nAuthority",
                                    "NPI (full)")),
         outcomes_s = case_when(outcomes == "covid_mask_own" ~ "(1) Wears mask",
                                outcomes == "covid_mask_others" ~ 
                                  "(2) Others should wear mask",
                                outcomes == "covid_mask_told" ~ 
                                  "(3) Told others to wear mask",
                                outcomes == "covid_vaccine_2" ~ 
                                  "(4) Get vaccinated"),
         is_sig = ifelse((estimate + 1.96*std.error) * (estimate - 1.96*std.error) > 0 ,
                         1, 0))

ggplot(res_weighted, aes(x = term_f, y = estimate,
             ymin = estimate - 1.96*std.error,
             ymax = estimate + 1.96*std.error, shape = wc,
             colour = factor(is_sig))) +
  geom_hline(yintercept = 0, colour = "grey80", alpha = 0.5, linetype = 2) +
  scale_colour_manual(values = c("grey80", "black"), guide = "none") +
  geom_linerange(position = position_dodge(width = 0.35)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  facet_wrap(~outcomes_s, ncol = 4) + 
  coord_flip() +
  labs(x = "", y = "Effect associated with 2SD difference (logit)",
       shape = "Model", caption = "Whiskers are 95% confidence intervals.\nStatistically significant (p < 0.05) coefficients in black.") +
  theme_minimal() +
  scale_shape_manual(values = c(16, 15)) +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", 
        strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom")

# uncomment to save
# ggsave(file = "weighted-model-coefs.pdf",
#        width = 12, height = 6,
#        device = cairo_pdf)

# summarize weighted vs not weighted comparison
res_unweighted$weights <- "No survey weights"
res_weighted$weights   <- "With survey weights"
res_all <- bind_rows(res_unweighted, 
                     res_weighted) |> 
  select(-model)

# comparative figure without controls (1)
ggplot(filter(res_all, wc == "without controls"), aes(x = term_f, y = estimate,
                         ymin = estimate - 1.96*std.error,
                         ymax = estimate + 1.96*std.error, shape = weights,
                         colour = factor(is_sig))) +
  geom_hline(yintercept = 0, colour = "grey80", alpha = 0.5, linetype = 2) +
  scale_colour_manual(values = c("grey80", "black"), guide = "none") +
  geom_linerange(position = position_dodge(width = 0.35)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  facet_wrap(~outcomes_s, ncol = 4) + 
  coord_flip() +
  labs(x = "", y = "Effect associated with 2SD difference (logit). No additional controls included.",
       shape = "Model", caption = "Whiskers are 95% confidence intervals.\nStatistically significant (p < 0.05) coefficients in black.") +
  theme_minimal() +
  scale_shape_manual(values = c(16, 15)) +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom")


# comparative figure with controls (2)
ggplot(filter(res_all, wc == "with controls"), aes(x = term_f, y = estimate,
                                                      ymin = estimate - 1.96*std.error,
                                                      ymax = estimate + 1.96*std.error, shape = weights,
                                                      colour = factor(is_sig))) +
  geom_hline(yintercept = 0, colour = "grey80", alpha = 0.5, linetype = 2) +
  scale_colour_manual(values = c("grey80", "black"), guide = "none") +
  geom_linerange(position = position_dodge(width = 0.35)) +
  geom_point(position = position_dodge(width = 0.35), size = 3) +
  facet_wrap(~outcomes_s, ncol = 4) + 
  coord_flip() +
  labs(x = "", y = "Effect associated with 2SD difference (logit). Controls included.",
       shape = "Model", caption = "Whiskers are 95% confidence intervals.\nStatistically significant (p < 0.05) coefficients in black.") +
  theme_minimal() +
  scale_shape_manual(values = c(16, 15)) +
  theme(text = element_text(colour = "grey20"),
        panel.grid.major = element_line(linetype = 3),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(hjust = 0),
        axis.text.x = element_text(colour = "grey20"),
        axis.text.y = element_text(colour = "grey20"),
        strip.text = element_text(colour = "grey20")) +
  theme(panel.background = element_rect(fill = "grey98", colour = NA)) +
  theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0),
        legend.position = "bottom")

# Session info for last run instance
# > sessionInfo()
# R version 4.4.0 (2024-04-24)
# Platform: aarch64-apple-darwin20
# Running under: macOS 15.1.1
# 
# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
# LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# time zone: Europe/Copenhagen
# tzcode source: internal
# 
# attached base packages:
#   [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] conflicted_1.2.0 survey_4.4-2     survival_3.5-8   Matrix_1.7-0     broom_1.0.6     
# [6] texreg_1.39.3    Hmisc_5.1-3      lavaan_0.6-19    psych_2.4.6.26   lubridate_1.9.3 
# [11] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2      readr_2.1.5     
# [16] tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1    tidyverse_2.0.0 
# 
# loaded via a namespace (and not attached):
#   [1] tidyselect_1.2.1  farver_2.1.2      fastmap_1.2.0     digest_0.6.35    
# [5] rpart_4.1.23      timechange_0.3.0  lifecycle_1.0.4   cluster_2.1.6    
# [9] magrittr_2.0.3    compiler_4.4.0    rlang_1.1.4       tools_4.4.0      
# [13] utf8_1.2.4        data.table_1.15.4 knitr_1.46        labeling_0.4.3   
# [17] htmlwidgets_1.6.4 bit_4.0.5         mnormt_2.1.1      withr_3.0.0      
# [21] foreign_0.8-86    nnet_7.3-19       stats4_4.4.0      fansi_1.0.6      
# [25] colorspace_2.1-0  scales_1.3.0      MASS_7.3-60.2     cli_3.6.2        
# [29] rmarkdown_2.27    crayon_1.5.2      ragg_1.3.2        generics_0.1.3   
# [33] rstudioapi_0.16.0 httr_1.4.7        tzdb_0.4.0        DBI_1.2.2        
# [37] cachem_1.1.0      splines_4.4.0     parallel_4.4.0    base64enc_0.1-3  
# [41] mitools_2.4       vctrs_0.6.5       hms_1.1.3         bit64_4.0.5      
# [45] Formula_1.2-5     htmlTable_2.4.2   systemfonts_1.1.0 glue_1.7.0       
# [49] stringi_1.8.4     gtable_0.3.5      quadprog_1.5-8    munsell_0.5.1    
# [53] pillar_1.9.0      htmltools_0.5.8.1 R6_2.5.1          textshaping_0.3.7
# [57] vroom_1.6.5       evaluate_0.23     pbivnorm_0.6.0    lattice_0.22-6   
# [61] backports_1.5.0   memoise_2.0.1     Rcpp_1.0.12       gridExtra_2.3    
# [65] nlme_3.1-164      checkmate_2.3.1   xfun_0.48         pkgconfig_2.0.3  